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ABSTRACT 

We study the gravitational collapse of a rotating supermassive star by means of a (3+1) hydrodynam- 
ical simulation in a post-Newtonian approximation of general relativity. This problem is particularly 
challenging because of the vast dynamical range in space which must be covered in the course of collapse. 
We evolve a uniformly rotating supermassive star from the onset of radial instability at R p /M = 411, 
where R p is the proper polar radius of the star and M is the total mass-energy, to the point at which the 
post-Newtonian approximation breaks down. We introduce a scale factor and a "comoving" coordinate 
to handle the large variation in radius during the collapse (8 < R p /Mo < 411, where Mq is the rest 
mass) and focus on the central core of the supermassive star. Since T/W, the ratio of the rotational 
kinetic energy to the gravitational binding energy, is nearly proportional to 1/R P for an n = 3 polytropic 
star throughout the collapse, the imploding star may ultimately exceed the critical value of T/W for 
dynamical instability to bar-mode formation. Analytic estimates suggest that this should occur near 
Rp/M ~ 12, at which point T/W ~ 0.27. However, for stars rotating uniformly at the onset of collapse, 
we do not find any unstable growth of bars prior to the termination of our simulation at R p /Mq ~ 8. 
We do find that the collapse is likely to form a supermassive black hole coherently, with almost all of 
the matter falling into the hole, leaving very little ejected matter to form a disk. In the absence of 
nonaxisymmetric bar formation, the collapse of a uniformly rotating supermassive star does not lead to 
appreciable quasi-periodic gravitational wave emission by the time our integrations terminate. However, 
the coherent nature of the implosion suggests that rotating supermassive star collapse will be a promising 
source of gravitational wave bursts. We also expect that, following black hole formation, long wavelength 
quasi-periodic waves will result from quasi-normal ringing. These waves may be detectable by the Laser 
Interferometer Space Antenna (LISA). 

Subject headings: Gravitation — gravitational waves — hydrodynamics — instabilities — relativity — 
stars: rotation 
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1. INTRODUCTION 

There is increasing evidence that supermassive black 
holes (SMBHs) exist at the center of all galaxies, and that 
they are the sources which power active galactic nuclei 
and quasars (Rees 1998; Macchetto 1999). For example, 
VLBI observations of the Keplerian disk around an object 
in NGC4258 indicate that the central object has a mass 
M ~ 3.6 x 1O 7 M and radius less than ~ 13pc. Also, large 
numbers of observations are provided by the Hubble space 
telescope suggesting that SMBHs exist in galaxies such as 
M31 (3 x 1O 7 M ), M87 (1 - 2 x 10 9 M Q ) and our own 
galaxy (2.5 x 1O 6 M ) (see for example, Macchetto 1999, 
for a brief overview). 

Although evidence of the existence of SMBHs is com- 
pelling, the actual formation process of these objects is 
still uncertain (Rees 2001). Several different scenarios have 
been proposed, some based on stellar dynamics, others on 
gas hydrodynamics, and still others which combine the 
processes. In one stellar dynamical scenario, a dense star 
system composed of compact stars becomes dynamically 
unstable to a collisionless, relativistic radial mode and un- 
dergoes catastrophic collapse to a SMBH (Zel'dovich & 
Podurets 1965; Shapiro & Teukolsky 1985a, b; Quinlan & 
Shapiro 1987). In an alternative scenario, massive stars 



build up within a dense cluster, following collisions and 
mergers of ordinary stars. After repeated collisions and 
mergers, supermassive stars (SMSs) are formed, and these 
become unstable to a hydrodynamic, relativistic radial 
mode (Iben 1963; Chandrasekhar 1964a, b; Feynman, un- 
published, as quoted in Fowler 1964) and eventually col- 
lapse to form SMBHs (Lee 1987; Quinlan & Shapiro 1990; 
Ebisuzaki et al. 2001). In still another scenario, massive 
halos of self-interacting dark matter in the early universe 
undergo the gravothermal catastrophe (secular core col- 
lapse), followed by catastrophic collapse to SMBHs once 
their cores become relativistic (Balberg, Shapiro, & Ina- 
gaki 2002; Balberg & Shapiro 2001). In a typical gas hy- 
drodynamical scenario, a contracting primordial gas cloud 
builds up sufficient radiation pressure to inhibit fragmen- 
tation, and the gas directly builds up a SMS (Sanders 
1970; Begelman & Rees 1978; Haehnelt, Natarajan, & Rees 
1998). Some mass inevitably will be shed but most of 
the matter is trapped during the ensuing gravitational col- 
lapse, forming a SMBH. At present, there is no definitive 
observation as yet which confirms or rules out any one of 
these scenarios. 

Here we focus on the collapse of a SMS. Baumgarte & 
Shapiro (1999) investigated the equilibrium, stability and 
quasi-static evolution of a SMS with uniform rotation. 



2 



SAIJO ET AL. 



During quasi-static evolution, uniform rotation can be 
maintained either by internal viscosity or magnetic brak- 
ing. They showed that the nondimensional ratios Rp/M , 
T/W and J/M 2 for all critical configurations at the on- 
set of collapse are universal numbers, independent of the 
history or mass of the star. Here R p is the proper polar 
radius, M is the gravitational mass (total mass-energy), 
T is the rotational kinetic energy, W is the gravitational 
binding energy and J is the angular momentum. They 
also pointed out the possibility of bar formation during 
catastrophic collapse prior to BH formation, assuming that 
the collapse is nearly homologous (T/W oc 1/R p ). They 
provided a crude analytic argument suggesting that the 
imploding star should pass the dynamical bar-mode in- 
stability point T/W ~ 0.27 at the radius R p /M ~ 12, 
but were not able to assess the consequences of this fact. 
New & Shapiro (2001) later investigated the quasi-static 
evolution of a SMS with differential rotation, assuming 
negligible viscosity and magnetic fields. They showed that 
in this case, bar formation prior to the onset of relativistic 
instability is inevitable. 

One of the primary observational missions for space- 
based detection of gravitational waves is the investigation 
of supermassive objects (Thorne 1998). Since the Laser 
Interferometer Space Antenna (LISA) will have long arms 
(10 6 km), the detector will be most sensitive in the low 
frequency band (10~ 4 ~ 10 _1 Hz). Potential sources of 
high signal to noise events in this frequency range include 
quasi-periodic waves arising from nonaxisymmctric bars 
in collapsing SMSs and the inspiral of binary SMBHs. In 
addition, the nonsphcrical collapse of rotating SMSs to 
SMBHs could be a significant source of burst and quasi- 
normal ringing radiation. In this paper we tract the col- 
lapse of a SMS by numerical simulation to investigate some 
of these possibilities. 

The growth of bars during the catastrophic collapse of 
stars has been studied previously for scenarios leading 
to core bounce, as in supernova core collapse. Rampp, 
Miiller, & Ruffert (1998) employed a two-component equa- 
tion of state with an effective adiabatic index larger than 
4/3, which causes bounce. They concluded that the star 
forms a bar in several dynamical timescales after it ex- 
ceeds the critical value for dynamical bar-mode instability 
(T/W > 0.27 for Newtonian incompressible equilibrium 
stars). But since the rotating core re-expands after the 
bounce and T/W falls to a lower value < 0.27, the non- 
axisymmetric instability does not have sufficient time to 
enter the nonlinear regime and the bar does not grow. 
Therefore, gravitational radiation is not generated effec- 
tively in this core bounce. Brown (2001) considered the 
same problem but chose a different equation of state to al- 
low the star more time (by a factor of 100) to reside in the 
unstable regime T/W > 0.27. Nevertheless, he too found 
that the dynamical instability was too weak to produce 
an appreciable bar and a large emission of gravitational 
waves. No one has considered the growth of bars during 
collapse without a bounce, which is the case for a SMS 
where the adiabatic index is close to 4/3. 

We take our initial stellar model to be a marginally un- 
stable SMS star near the critical point, R p /M ~ 430. We 
treat the gas adiabatically, since for sufficiently massive 
stars neither photon nor neutrino losses are dynamically 



significant (Linke et al. 2001). We take the adiabatic in- 
dex to be 4/3, appropriate for a radiation-pressure dom- 
inated SMS, and construct a critical, uniformly rotating 
polytrope with index n w 3 for our starting point. Our 
goal is to determine the final outcome of the collapse. We 
want to address the following questions: Does a SMBH 
definitely form following the catastrophic collapse? Is the 
collapse coherent or does the central region collapse first, 
followed by the gradual accretion of the envelope? Does 
the collapsing configuration fragment? Does a disk form? 
Does a rotating bar form during the collapse? 

We use a post-Newtonian (PN) hybrid hydrodynami- 
cal code in (3+1) dimensions to tract SMS collapse. Our 
adopted hybrid scheme is relativistically exact for static 
spherical spacetimes (Shibata, Baumgarte, & Shapiro 
1998). The onset of radial instability occurs when T/W -C 
0.1, so our initial equilibrium spacetime is very nearly 
spherical. Locating the onset of radial instability in a 
SMS requires the presence of nonlinear gravitation to at 
least 2PN order (Zel'dovich & Novikov 1967; Baumgarte 
& Shapiro 1999). For these reasons, the nonlinearity cap- 
tured in our hybrid scheme, which extends beyond 1PN, 
is essential to treat this problem. Of course, it is neces- 
sary to use a fully general relativistic code to follow the 
final implosion of the matter into a black hole (BH) and to 
reliably determine the gravitational waveforms. (Shapiro 
& Teukolsky (1979) followed SMS collapse in full general 
relativity for a nonrotating configuration with a spheri- 
cal [1+1] code). However, a fully relativistic (3+1) code 
capable of handling the large dynamic range spanned by 
SMS collapse is not yet available. Fortunately, since our 
initial configuration is nearly Newtonian, we can use our 
hybrid scheme to track most of the implosion up to the 
point where the formation of a BH is likely. Our hybrid 
scheme is also adequate to address most of the questions 
raised above, at least in a preliminary fashion. 

This paper is organized as follows. In Sec. 2 we present 
the basic equations of our PN formulation in "comoving" 
coordinates. We demonstrate the ability of our code to dis- 
tinguish stable from unstable stars in Sec. 3. We discuss 
our numerical results for catastrophic collapse in Sec. 4. In 
Sec. 5 we summarize our findings. Throughout this paper, 
we use geometrized units (G = c = 1) and adopt Cartesian 
coordinates (x, y, z) with the coordinate time t. Greek and 
Latin indices run over (t, x, y, z) and (a;, y, z), respectively. 

2. NUMERICAL METHOD AND KEY EQUATIONS 

In this section, we briefly derive the (3+1) hybrid PN 
relativistic hydrodynamic equations (Shibata, Oohara, & 
Nakamura 1997; Shibata, Baumgarte, & Shapiro 1998) in 
"comoving" coordinates. We solve the fully relativistic 
equations for hydrodynamics, but neglect some higher- 
order dynamical PN terms in the Einstein field equations. 
Note that this approximation gives the exact solution for 
a static spherical spacetime. To track the collapse over the 
vast dynamic range from > 410M down to a few M and 
to investigate the central core at late times, we require 
a suitable comoving coordinate system. Such a coordi- 
nate choice is possible because in Newtonian gravity, an 
n = 3 spherical polytrope collapses homologously (Goldre- 
ich & Weber 1980). This special behavior does not strictly 
hold in general relativity (Shapiro & Teukolsky 1979), nor 
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does it apply to a rotating configuration, but it holds ap- 
proximately for much of the implosion of a critical SMS 
configuration, since it is nearly Newtonian and slowly ro- 
tating. Therefore, we can construct a "comoving" frame, 
subtracting the mean "Hubble" flow from the local ve- 
locity, to follow most of the collapse with sufficient grid 
resolution. 

First we construct a PN metric in the "comoving" frame. 
The line element in the PN formalism is written as (Chan- 
drasekhar 1965; Blanchet, Damour, & Schafer 1989) 

ds 2 — g^dx^dx" 



where 



= (- 



■a 



■ I3 k /3 k )dt 2 + 2f3 k dx k dt + ^8 i3 dx l dx\ (1) 



where a, p, and ip is a lapse function, shift vector, and 
conformal factor, respectively. Therefore, as in the study 
of BH formation in the Friedman Universe (Shibata & 
Sasaki 1999), we define the "comoving" frame by 

x l = dx\ (2) 

where x % is the spatial coordinate in the "comoving" frame, 
and a is a scale factor which only depends on t. Hereafter, 
A will represent the quantity A as measured in "comov- 
ing" coordinates. Note that we do not change the time 
coordinate. We write the PN line element in "comoving" 
coordinates as 

ds 2 — g^dx^dx" 

= (-a 2 + j3 k (3 k )dt 2 + 2$ k dx k dt + d 2 tp 4 5 i3 dx l dx j (3) 

By matching the two line elements (eqs. [1] and [3]) and 
using equation (2), we identify 

a = a, (4) 

I3 k 



P 



Hot* 



(5) 

ip = ip, (6) 

where H = a/a. Note that, among the geometric quanti- 
ties, we only need to adjust the shift vector. 

For a perfect fluid, the energy-momentum tensor in "co- 
moving" coordinates is written as 

f^v = P (l + e + UpU v + Pg^, (7) 

where p, e, and P are rest mass density, specific internal 
energy density, and pressure, respectively. We adopt a 
T-law equation in the form P = (T — l)pe in this paper. 

In the "comoving" frame, the continuity equation, en- 
ergy equation, and Euler equation with T-law equation of 
state including artificial viscosity are written as 



— (a J p*^) + ^jj( a P*UiV J ) 

d ,„ „ s ,o ~t da 



(8) 



(9) 



- dp 



'^^W^ + Pvis) ~ a " p * a " dx 1 + a0p * flj ' W 



(10) 



v = — > 



u = (1 + Te)u t , 
iii = (1 + Te)ui. 



(11) 
(12) 

(13) 
(14) 
(15) 



Following Shibata (1999) we include artificial viscosity in 
our system as 



3 r 

'vis — \ (era ip ) 



(16) 



0, 



for Sv > 0, 



with Sv = 25xdiV l 1 C v is — 1, and Sx is the step size of the 
"comoving" grid. 

Gravitational field equations in the PN approximation 
are derived from the Hamiltonian constraint, momentum 
constraint and the maximal time-slicing condition (Shi- 
bata, Oohara, & Nakamura 1997; Shibata, Baumgarte, & 
Shapiro 1998). The equations in the comoving frame are 
written as 



A-0 = -27ra 2 -0 5 /5 H , 

A(cn/0 = 27raa 2 ?A 5 (pH + 25), 



SijAp + -didjP = 16naJi, 



(17) 
(18) 

(19) 



where p H = n^T"", J t = -n^T"", S = h^h lv T^ , 

= (—a, 0,0,0), h^u = g^ u + n^n v , and A is the flat 
Laplacian measured in the "comoving" frame. 

We use the asymptotic fall-off behavior for metric quan- 
tities at large radius in order to set an appropriate bound- 
ary condition at the grid edge (Shibata et al. 1998). Only 
the rescaled shift vector p requires special consideration. 
Although P has a simple fall-off at the large radii, falling 
like <~ 0(r~ 2 ), P diverges like ~ Hx l . In order to elim- 
inate the divergent behavior of p, we introduce a new 
vector W l as 



W i = P - Hx\ 



(20) 



In terms of W % , the momentum constraint equation (eq. 
[19]) can be rewritten as 



5 l3 AW j + \didjW j = 16naJi. 



(21) 



From a computational point of view, it is useful to split 
this momentum equation into one vector Poisson equation 
and one scalar Poisson equation using the technique of 
Shibata (1997): 

SijWi = 4B, - ±[da + di(B k x k )], (22) 

where E>i and x satisfy 

AB i = 4naJ i , (23) 
Ax - -AitaJiX 1 . (24) 
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We also need to rewrite the Euler equation [eq. (10)] using 
W* as 

— (a 2 p*Ui) + j—^p^u.v ) 



d 



-a z a^—AP + P wis ) 
ox 1 

2 1 + Te . 2 dip 

+ 2 « P* -rrr- [(««*) -i]^' 



dx' 1 



- 



ox 1 



where 



(25) 



(26) 



The fluid velocities in the two frames are related through 
equation (2) by 

v l = a(v i +Hx l ). (27) 

Finally, we must decide how to choose H. As long as 
homology holds during the collapse, we can always sub- 
tract the "Hubble" flow by choosing the appropriate H 
(v r = 0). This is why the star docs not contract during 
the collapse in the "comoving" frame. Since rotating col- 
lapse does not strictly hold homology due to the effects 
of rotation and PN gravity, success in exploring the late 
collapse depends on the location where we choose to sub- 
tract the "Hubble" flow. Though ann = 3 polytropic star 
has a large envelope, most of the mass is located in the 
central core. Therefore we set H around the mean radius 
according to 



H 



-W x 



where the mean radius is defined as 



/x x 



7 m 



' / p*x 2 d 3 x 



(28) 



(29) 



We experiment by using different radii in the neighbor- 
hood of r m for measuring H and choose the one which 
yields the most accurate evolution. Note that for homolo- 
gous Newtonian collapse of an n — 3 spherical polytrope, 
H defined above is independent of the radius at which it 
is evaluated. 

We monitor several global quantities of the system dur- 
ing the collapse. The gravitational mass M, rest mass Mo, 
proper mass M p , angular momentum J, rotational kinetic 
energy T, and gravitational binding energy W of the ro- 
tating star are defined as 

1 ' VtydSi 



M = 



2ir 



- J a 3 (p + pe + P)(mt*) 2 -P 
J pu t \f^gd?x — J a 3 p r d 3 x, 
J pu t (l + e) v ^d 3 x 
J a 3 p*(l + e)d 3 x, 
l - ( (xK* - yK^dSi 

>7T J r — >oq 



ip 5 d 3 x, 



8tt 



d 3 (xJ v - yj x )tp 6 d 3 x 



(30) 
(31) 

(32) 

(33) 



1 



W=\M V + T 



2 
M\, 



QdJ = - arQ(xJ y 



yJ x )^d 3 x, (34) 
(35) 



where fi is the angular velocity of the star, 



yv 



x 2 + y 2 



(36) 



Note that in PN gravity, M and J are not strictly con- 
served, especially when the field becomes strong and 
higher order corrections to the gravitational field become 
important. However, these quantities should be well con- 
served during those early and intermediate epochs when 
gravity remains weak. Also J defined above should be 
strictly conserved whenever the system is axisymmetric. 
Since our system is nearly axisymmetry, prior to any bar 
formation, monitoring J conservation enables us to check 
the accuracy of our numerical results. We also define the 
quadrupole moment hj as 



and the nonaxisymmetric distortion parameter r\ as 

1 xx lyy 



^xx H~ 



(37) 



(38) 



As we choose a polytropic equation of state for our initial 
data, with P = np 1+1 ^ n , where n w 3 and n is a constant 
appropriate for thermal radiation pressure of constant spe- 
cific entropy (see, e.g., Baumgarte & Shapiro 1999, eq. 
[3]). For numerical purposes it is convenient to rescale all 
quantities with respect to n. Since n 11 / 2 has dimensions 
of length, we introduce the following nondimensional vari- 
ables 



M_=K~ n / 2 M, R = n- n / 2 R, J = K~ n J, 

n = K n / 2 n, p = K n P . 



(39) 



Therefore our numerical results apply to arbitrary mass, 
where we simply need to adjust k according to 



M 



10 6 Mr 



,n/2 



1.48 x 10 6 km 



(40) 



to apply our result to arbitrary mass. Henceforth, we 
adopt nondimensional quantities, but omit the bars for 
convenience. Equivalently, we set k = 1. 



3. STABILITY ANALYSIS 

Before simulating rotating SMS collapse in our "comov- 
ing" scheme, we first assess the stability and accuracy of 
our (3+1) PN hybrid code. For this purpose, we study the 
stability of stars along an equilibrium sequence of fixed en- 
tropy (fixed k) using our dynamical code and compare the 
results with the exact turning point criterion. We recall 
that the onset of radial instability occurs at the turning 
point in a plot of M vs. p c along the sequence. 
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3.1. Spherical Symmetric Case 

In general relativity and PN gravitation, all nonrotating 
spherical polytropcs with index n — 3 are radially unstable 
to collapse and have no turning points along their equilib- 
rium sequences. To test the ability of our code to distin- 
guish stable from unstable spherical stars, we must study 
polytropes with a slightly smaller n; we set n = 2.96 for 
this purpose, for which the equilibrium structure remains 
close to that of a radiation-pressure supported SMS. We 
perturb the initial equilibrium at t = by slightly decreas- 
ing the pressure and following the subsequent evolution. 
More specifically, we decrease the pressure constant k ac- 
cording to k — ► 0.99k. We use modest grid sizes (typically 
99 x 99 x 50) 1 for this test, covering the star with 81 grid 
points across its diameter. 
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Fig. 1. — Probing the dynamical stability of a spherical SMS 
with n = 2.96. Here, p c is the central density of the equilibrium 
spherical star. Filled circles and crosses represent unstable stars, 
while open circles represent stable stars according to our dynamical 
calculation. A cross indicates that the star is actually stable ana- 
lytically according to the turning point criterion. The radii of the 
8 marked stars are R/M = 150, 200, 300, 750, 1007, 1499, 2001, 
4003, where the sequence starts at the right side of the figure at 
the highest central densities. With the adopted grid resolution, our 
code can distinguish stable from unstable stars to within 1% of the 
maximum gravitational mass. 

Figure 1 is a summary of our stability code calibration 
for a spherical SMS. We find that with the above grid reso- 
lution, our code can distinguish stable from unstable stars 
to within 1% of the maximum gravitational mass. Figure 
2 compares the dynamical behavior of stable and unstable 
stars. Note that we use the unit of time as tp = J R^/M 
where R c is the initial equatorial radius (R c = R for a 
spherical case). For dynamically unstable stars the cen- 
tral density increases monotonically, while for stable stars 
it oscillates. 

3.2. Rotating Case 

To assess the ability of our code to distinguish stable 
from unstable stars with rotation, we consider an equilib- 
rium sequence of uniformly rotating stars of fixed angular 
momentum J (J/M 2 = 0.644 at the turning point). While 
the turning point criterion strictly identifies the onset of 
secular instability (Friedman, Ipser, & Sorkin 1988), 

1 Our code has equatorial plane symmetry. 




2 4 6 



t/t D 

Fig. 2. — Evolution of the central densities of the stars plotted 
in Fig. 1. Here trj = Rc\> Rc/M, where Re is the initial equatorial 
radius. Curves are drawn for stars which are unstable both numer- 
ically and according to the turning point criterion (solid), unsta- 
ble numerically but stable according to the turning point criterion 
(dashed), and stable both numerically and according to the turning 
point criterion (dash-dotted). 

the point of onset of dynamical instability nearly coin- 
cides with the secular instability point (Shibata, Baum- 
garte, & Shapiro 2000). We adopt the same polytropic 
index and grid resolution as in the spherical simulations 
reported above. We decrease the initial pressure as in the 
spherical case (k — > 0.99k). 
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Fig. 3. — Probing the dynamical stability of a rotating SMS with 
n = 2.96, J = 10. Circles and crosses have the same meanings as 
in Fig. 1. The radii of the 5 marked stars are R/M = 252, 285, 
419, 537, and 785, where the sequence starts at the right side of the 
figure at the highest central densities. Note that the solid line shows 
a constant J sequence with J = 10, while the dashed line represents 
the spherical equilibrium sequence (Fig. 1). With the adopted grid 
resolution, our code can distinguish stable from unstable stars to 
within 1% of the maximum gravitational mass. 

Figure 3 summarizes our dynamical stability analysis for 
the rotating SMS. We conclude that with the adopted grid 
resolution, our code can distinguish stable from unstable 
rotating stars to within 1% of the maximum gravitational 
mass. Figure 4 shows the evolution of the central density 
for stable and unstable rotating stars. 
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Fig. 4. — Evolution of the central densities of the stars plot- 
ted in Fig. 3. Curves are drawn for stars which are unstable both 
numerically and according to the turning point criterion (solid), un- 
stable numerically but stable according to the turning point criterion 
(dashed), and stable both numerically and according to the turning 
point criterion (dash-dotted). 

4. CATASTROPHIC COLLAPSE 

In this section, we simulate catastrophic collapse, focus- 
ing on the possibility of bar formation and on the final fate 
of the central object. First, we examine spherical SMS col- 
lapse in PN hybrid gravitation, to verify numerically that 
we can recover homology over the expected range of radii 
and that we can resolve the late stages of the collapse as 
long as general relativistic gravity is not too strong. Next, 
we examine the collapse of a rapidly rotating n = 1 poly- 
trope to check that our code can reproduce bar formation 
when the collapse is followed by a bounce. Finally, we in- 
vestigate the collapse of a uniformly rotating SMS, which 
is the main result of this paper. We discuss the final fate 
of the SMS collapse, the formation of a SMBH and the 
possibility of circumstellar disk formation and mass loss. 

4.1. Spherical Collapse 

We examine spherical SMS collapse to check the range of 
radii over which homology holds during the collapse in PN 
hybrid gravitation. For the collapse of an n = 3 spherical 
Newtonian star, self-similarity is strictly preserved (Gol- 
dreich & Weber 1980). In full general relativity the col- 
lapse of an n — 3 spherical star maintains self-similarity 
down to the radius R/M < 80 (Shapiro & Teukolsky 
1979). Therefore, we expect that there is some region over 
which self-similarity is maintained in PN hybrid gravita- 
tion. Because n = 3 spherical stars are marginally un- 
bound, the system may expand rather than collapse for 
this polytropic index when we integrate with finite grid 
resolution. To guarantee collapse, we choose a slightly 
lower index, n — 2.96 for our computation. Since this com- 
putation is only to reconfirm homologous behavior above 
a minimum stellar radius and to probe how our code sig- 
nals the final fate of spherical collapse (a BH), we choose 
a modest initial radius R/M = 150. By performing this 
spherical collapse in 3D instead of ID, this simulation does 
provide a useful testbed for our subsequent rotating col- 
lapse simulation. The parameters characterizing our initial 
equilibrium star are summarized in Table 1. We acknowl- 
edge that our hybrid PN code cannot accurately handle 
BH formation, which is a fully general relativistic 



Table 1 
Parameters for the initial 
spherical equilibrium sms. 



Pc 


M 


R/M 


2.56 x 10- y 


3.80 


150 



phenomenon. However, we expect to be able to follow 
collapse far enough to ascertain whether or not BH for- 
mation is the likely outcome. 

We plot density profiles for the collapsing star at selected 
times in Figure 5. We decrease the initial pressure of the 
equilibrium star by 1% (k — > 0.99k) to induce collapse. 
We see that the collapse is nearly homologous up to a time 
t/tn ~ 2.5. The radius at that time is R/M ~ 80. Note 
that the central density exceeds 20 times its initial value 
at this point. Our result is consistent with that of Shapiro 
& Teukolsky (1979), who found that homology was main- 
tained until R/M <~ 80. We do not expect homologous 
evolution beyond this point, due to deviations from New- 
tonian behavior. We continue the collapse to t = 2.70tD, 
by which time the lapse has decreased to a c ~ 0.2 and the 
mean radius is about r m /M ~ 2.0 (recall that the horizon 
radius for a static BH in the adopted isotropic coordinate 
system is r/M — 0.5). Given the plunge in the lapse, we 
may safely guess that a BH will form after further collapse. 



1 



0.8 




Fig. 5. — Density profiles at selected times during spherical SMS 
collapse. Solid, dotted, dashed, long-dashed lines show profiles at 

t/t D = 0, 2.04, 2.52, 2.68 and (p c /pf \ R/M ) = (1, 146) , (3.05, 
131), (16.5, 76.0), (770, 24.8), respectively. Note that when a profile 
overlaps the solid line, the system is strictly homologous. Homol- 
ogy is maintained in this spherical collapse up to the time at which 
R/M ~ 80. 

4.2. Core Collapse and Bounce 

One scenario that leads to bar formation is core bounce 
following the collapse of a rotating star. Simulating this 
phenomenon tests the ability of our code to detect a bar- 
mode instability. Two previous numerical studies focus on 
nonaxisymmetric instabilities arising in the core-collapse 
supernova (Rampp et al. 1998; Brown 2001). The effec- 
tive r in both cases is larger than 4/3 at the time of max- 
imum compression, which causes the bounce. They con- 
clude that once the star exceeds the value T/W ~ 0.27, bar 
formation takes place after several dynamical timescales. 
We use a similar scenario to test our code by reproducing 
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this result. We adopt a stiff n = 1 (r = 2) polytropic 
model to guarantee a core bounce following a modest de- 
crease in initial radius. Since the change in radius is not 
large, we do not need to install a scale factor in this cal- 
culation. As core bounce will generate a shock between 
the infalling and outgoing matter, we must utilize artifi- 
cial viscosity in our code. 

The initial data for this case is summarized in Table 2. 
We choose initial data near the mass-shedding sequence 
in order to make the initial value T/W as large as possi- 
ble. We also choose the initial equilibrium star so that the 
equatorial radius remains in the PN regime, R/M > 20. In 
order to achieve an appreciable implosion prior to bounce 
we drastically deplete the initial pressure by decreasing n 
according to 

k->0.09xk. (41) 

We also impose a slight triaxial density perturbation on 
the equilibrium star to induce m — 2 bar formation 
(Rampp et al. 1998), 



^(equilibrium) 



l 



y 



(42) 



where S = 0.1. We employ a grid size (239 x 239 x 120) 
which covers the initial equilibrium star with 161 points 
across the equatorial diameter. In fact, the proper equa- 
torial radius varies between 30 < R e /M < 104 during the 




collapse, so the resolution of the central region is always 
adequate during the evolution. 



0.3 

0.2 

0.1 


0.1 
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Fig. 6. — Growth of the bar in rotating core collapse and bounce. 
The quantity T/W exceeds the critical dynamical instability point 
~ 0.27 during the collapse and settles down below this value to 
T/W ~ 0.19. The deformation parameter r\ grows exponentially 
during the early evolution following maximum compression and os- 
cillates about a finite value at late times. Note that the behavior of 
T/W is quite similar to Fig. 3 of Rampp, Miiller, & Ruffert (1998). 
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Fig. 7. — Density contours p* in the equatorial plane at selected times during the collapse and bounce of a rotating star. Snapshots are 
plotted at (t/tu, Pt) = ( a ) ( 4 - 90 x 10_5 > 5 - 353 x iO" 3 ), (b) (1-41, 1.410 x lO" 1 ), (c) (3.43, 1.435 x 10" 1 ), (d) (4.70, 1.478 x lO" 1 ). The 
contour lines denote densities p* = p* X i/16 (* = !,•••, 15). 
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Table 2 

Parameters for the initial rotating equilibrium n = 1 polytrope. 



p c M J R p /M 


J/A-P 


T/W 


Rp/ R e 


5.00 x lQ- :i 1.34 x 10~ a 4.43 x 10~ 4 70.4 


1.91 


8.18 x 10~ 2 


0.675 



The bar diagnostics are shown in Figure 6. Our results 
are similar to those of Rampp et al. (1998), especially the 
behavior of T/W (Fig. 3 of Rampp et al. 1998). Figure 
6 shows the evolution of the nonaxisymmetric distortion 
function. Its growth signifies bar formation. 




-30 -20 -10 10 20 30 

x/M fl 



Fig. 8. — Density contours p* in the meridional plane at selected 
times during the collapse and bounce of a rotating star. The times 
and contour levels are the same as in Fig. 7. 

We plot density contours at selected times in Figures 
7 (equatorial plane) and 8 (meridional plane). A modest 
spiral arm forms within a dynamical time after the star ex- 
ceeds T/W = 0.27 (Fig. 7 [b]). After that, the spiral arm 
is twisted around the star and settles down to a weakened 
bar structure at the central core. This whole picture 



is qualitatively similar to the result reported in Rampp et 
al. (1998). 

From the density snapshots and the behavior of distor- 
tion function, we conclude that our code has successfully 
reproduced (weak) bar formation during core collapse and 
bounce. We thereby confirm that the code is capable of 
identifying bars when the physical conditions make it pos- 
sible for them to form. 

4.3. Rotating SMS Collapse 

Consider an overview of SMS evolution (Baumgarte & 
Shapiro 1999). Cooling and contraction of a rotating SMS 
will ultimately spin it up to the mass-shedding limit. After 
that, the SMS contracts secularly along the mass-shedding 
sequence as it cools, slowly losing mass and maintaining 
uniform rotation via viscosity and/or magnetic braking 
(Zel'dovich & Novikov 1967; Shapiro 2000). Upon reach- 
ing the onset of radial instability, the star will collapse 
catastrophically and form a BH, or a flattened rotating 
disk, or some combination thereof. It is this catastrophic 
collapse which we wish to follow with our dynamical code. 

We summarize the parameters of our initial uniformly 
rotating star in Table 3. We slightly perturb the initial 
equilibrium state according to 

k -> 0.99k, (43) 
P = p (°q uilib " um ) + S x2 ~v 2 ^ , ( 44 ) 

where 6 = 0.1. We slightly decrease k in order to de- 
plete the pressure and initiate the collapse. We install a 
triaxial density perturbation to provide the seed for bar 
formation, if the physical situation should lead to unsta- 
ble growth. We adopt a grid size (239 x 239 x 120), so 
that the star is initially covered by 161 points across the 
equatorial diameter. We evolve the rotating SMS up to 
the point at which the PN approximation breaks down. 

We show the evolution of the central lapse function in 
Figure 9. The figures shows that we can follow the col- 
lapse from the Newtonian regime where a c ~ 0.99 to the 
relativistic regime where a c ~ 0.3. The rapid plummet of 
a c below 0.3 at late times indicates that a BH will form 
as an immediate consequence of collapse. 



Table 3 

Parameters for the initial rotating equilibrium SMS. 





Pc 


M 


J 


Rp/M 


J/M' z 


T/W 


R p /R e 


Our initial data 


8.00 x 10" y 


4.57 


20.0 


411 


0.960 


8.85 x 10-* 


0.675 


Critical Value a 


7 x 10~ 9 


4.57 


20.3 


427 


0.97 


8.99 x 10~ 3 


0.664 



a Baumgarte & Shapiro (1999) 
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We plot the mass and angular momentum during the 
evolution in Figure 10. Our hybrid PN scheme conserves 
Mo, provided that the adopted cutoff density in the am- 
bient atmosphere is negligible 2 and that there is no mass 
loss from the grid. Our fully relativistic expressions for M 
and J, on the other hand, are not necessarily conserved in 
our hybrid PN scheme. We therefore chose to normalize all 
our length and time units in terms of M . Toward the end 
of our simulation, a small amount (< 4%) of the matter 
leaves the computational grid, leading to the loss of rest- 
mass as shown in Figure 10. This figure also demonstrates 
that M and J are conserved to similar levels, indicating 
that relativistic effects, which might cause larger devia- 
tions, are small. We show the evolution of the scale factor 
in Figure 11. This plot indicates that physical grid that 
we cover shrinks during the collapse from x max /M ~ 900 
(a = 1) to x max /Mo ~ 50 (a = 0.06) where x max is the 
edge of our numerical grid. 



We monitor the bar-mode diagnostics in Figures 12 and 
13. The amplitude of the deformation function decreases 
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Fig. 9. — Evolution of the central lapse function during rotating 
SMS collapse. 
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Fig. 10. — Evolution of rest mass Mo, gravitational mass M, 
and total angular momentum J during rotating SMS collapse. Note 
that Mo should be strictly conserved during the evolution within 
a numerical error and it is indeed conserved within 4% error (see 
text). 



Fig. 11. — Evolution of the scale factor during rotating SMS 
collapse. The collapse of the scale factor indicates that the physical 
size of our grid contracts from R/Mq ~ 900 (d = 1) to R/Mq ~ 50 
(a = 0.06). 
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Fig. 12. — Evolution of bar-mode diagnostics during rotating 
SMS collapse. The deformation parameter rj does not grow expo- 
nentially even in the final stage. 
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2 In order to regularize the fluid equations near the surface of the 
star it is convenient to put an atmosphere in the vacuum. We do 
this by setting a cutoff density below which p* may not drop. We 
define the surface of the star to be the place where p* = p* = 10p* ut , 
where p* ut /p* « 10" 9 . 



Fig. 13. — Growth of T/W during collapse. The solid line shows 
the results of our simulation (the same T/W as in Fig. 12). The 
dash line assumes T/W oc l/r m during the collapse; this scaling 
would hold if the collapse were spherical. Note that r m (0) is the 
mean radius at t = 0. 
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during the evolution, hence we conclude that a nonaxisym- 
metric bar does not grow prior to BH formation. Even 
though T/W may exceed « 0.27 in the final stages, there 
is not sufficient time for growth prior to the appearance of 
a hole. The reason for the absence of an instability is that 
T/W does not grow as rapidly as r" 1 during the late stage 
of collapse, as it would remain spherical (see Fig. 13). 

We show the density profile in the equatorial plane in 
Figure 14 and in the meridional plane in Figure 15. The 
ability of our scale factor implementation to resolve the 
matter distribution even as it becomes increasingly com- 
pact during the implosion is evident from these snapshots. 
We find no indication of the formation of a circumstellar 
disk with significant mass by the termination of our sim- 
ulation. In fact, the fraction of the rest mass outside a 
sphere of radius r/Mo — 7.0 is 26% and outside the sphere 
of r/Mo = 28.0 is 10%. Accordingly, most of the mass is 
concentrated in the center and is collapsing inward when 
we terminate our integration. Note that by employing a 
density cutoff, we are not reliably resolving the very out- 
ermost region. But we note that even with a cutoff, Mo is 
conserved to 96% accuracy (Fig. 10). We also note that 
the overall picture is not affected by a change of the cutoff 
value or extension of the grid size. We thus conclude that 
the rotation cannot provide sufficient centrifugal support 
in the bulk of the envelope to counter gravity and form a 
significant disk. 



Though our computation is terminated when the lapse 
drops below a c ~ 0.3, we can still infer the final fate of the 
collapse from examination of the velocity profile of the star 
(Figs. 16 and 17). The growth of an appreciable inward ra- 
dial component of the velocity field strongly suggests that 
immediately after the time we terminate the integrations 
the bulk of the matter will cross the event horizon of the 
nascent BH in a dynamical timescale as measured at the 
center of the star. 

Though the newly formed BH acquires the bulk of 
the mass in a coherent implosion, it does not obtain all 
of the mass and angular momentum. Lingering gaseous 
fragments in the outermost envelope containing negligible 
mass but nonnegligible angular momentum are not fol- 
lowed in our simulation, which focuses on the imploding 
massive bulk of the star. These fragments may accrete on 
a longer timescale (or even escape), but we cannot track 
their evolution with our current calculation. 

Figure 18 shows the angular momentum distribution 
during the evolution. During the collapse, the central lay- 
ers begin rotating faster than the surface layers and the 
configuration, uniformly rotating initially, acquires appre- 
ciable differential rotation. 

If all of the initial mass-energy and angular momentum 
are consumed by the final BH, it will be rapidly rotating 
with a/M « 1 (see Table 3). Although we cannot follow 
the final formation and growth of the BH, our PN simula- 
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Fig. 14. — Density contours p t in the equatorial plane at selected times during rotating SMS collapse. Snapshots are plotted at (t/tjy, p*, 
d) = (a) (5.0628 x 10" 4 , 8.254 x 10~ 9 , 10~ 7 ), (b) (2.50259, 1.225 x 10~ 4 , 10~ 5 ), (c) (2.05360, 8.328 x 10" 3 , 5.585 x 10~ 7 ), (d) (2.50405, 
3.425 X 10~ 2 , 1.357 X 10 -7 ), respectively. The contour lines denote densities p* = p* X c/f 1-1 / 1 ^) (i = 1, • • ■ , 15). 
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tions suggest that the final a/M may be slightly lower, 
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due to the loss of angular momentum carried by gas orbit- 
ing near the equator. 

5. DISCUSSION 

We follow rotating SMS collapse from the onset of radial 
collapse at R p /Mq ~ 411 to the point where the PN ap- 
proximation breaks down (R p /Mq ~ 8). The challenge of 
covering this large a dynamic range is met by introducing 
a scale factor and a "comoving" coordinate system which 
takes advantage of the homologous nature of the initial 
collapse. 

Collapse of a uniformly rotating, relativistically unsta- 
ble SMS is coherent and leads to the formation of a SMBH 
containing the bulk of the mass of the progenitor star. It 
is interesting to contrast our results with those of Loeb & 
Rasio (1994), who treat the isothermal (r = 1) collapse 
of initially homogeneous, uniformly rotating, low entropy 
clouds via smooth particle hydrodynamics (SPH) simula- 
tions. They find considerable fragmentation into dense 
clumps, and disk formation containing ~ 5% of the mass. 
They conclude that a seed BH will form at the center and 
that it likely will grow gradually by accretion. 



Fig. 15. — Density contours p* in the meridional plane at selected 
times during rotating SMS collapse. The times, the central densities 
and contour levels are the same as in Fig. 14. 
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Fig. 16. — Velocity field in the equatorial plane at selected times during rotating SMS collapse. The time for each snapshot is the same 
as in Fig. 14. Note that the velocity field is drawn in physical coordinates and the velocity arrows are normalized as indicated in the upper 
right hand corner of each snapshot. 



12 



SAIJO ET AL. 



We find no evidence of bars prior to BH formation, so 
that the collapse is largely axisymmetric. As a result, little 
angular momentum can be radiated away by gravitational 
waves. 

From the coherent, axisymmetric nature of the implo- 
sion we conclude that the collapse of a SMS, rotating uni- 
formly at the onset of collapse, is a promising source of 
gravitational wave bursts. We can estimate the strength 
and the frequency of the wave burst emitted from this 
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Fig. 17. — Velocity fields in the meridional plane at selected 
times during rotating SMS collapse. The time for each snapshot 
is the same to Fig. 14. We omit Fig. 17 (a) because there is no 
velocity in the meridional plane at t = 0. 



rotating collapse. The characteristic burst frequency is 
given by the dynamical timescale of the star at the time 
of BH formation, 
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The wave amplitude can be estimated by employing the 
Newtonian quadrupole formula according to 



^burst ^ i 



MR 2 f 2 



burst 



2 x 10 



-19 f M 



M 1 M 

where Q is the quadrupole moment of the star and d is the 
distance from the observer. We set R/M — 5, a charac- 
teristic mean radius during BH formation. Since the main 
targets of LISA are gravitational radiation sources between 
10 -4 and 10 _1 Hz, it is possible that LISA can search for 
the burst waves accompanying rotating SMS collapse and 
formation of a SMBH. 

In the absence of bar formation, SMS collapse will not 
produce quasi-periodic waves prior to SMBH formation. 
However, such waves will be generated by the nascent 
BH via quasi-normal mode ringing. The characteristic 
frequency /qnm and strength /iqnm of this radiation in 
rotating star collapse are (Thorne 1987; Shibata, Shapiro, 
& Uryu 2001) 



9x 10" 



0.02 



(a) 



a 



6x10 - 



3 x 10" 



a 




-500 -250 
0.1 



250 




500 -40 
-n 0.15 



-20 



20 



40 



(d) 



0.1 - 



0.05 




^ 1 
-30 -20 -10 10 20 30 

x/M 



Fig. 18. — Angular momentum distribution along the x-axis at selected times during rotating SMS collapse. The time for each snapshot 
is the same as Fig. 14. 
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where A£?gw is total radiated energy. Here we use the 
I = m = 2 quasi-normal mode frequency wqnm for a Kerr 
BH with a/M = 0.9, for which Mwqnm = 0.7 (Leaver 
1985). The efficiency of the radiated energy could be less 
than AE GW /M < 7 x 10~ 4 for a rotating collapse (Stark 
& Piran 1985). From the above estimate, gravitational 
waves from the vibration of a newly formed SMBH reside 
within the sensitivity limits of LISA. 

Since we assume that the equilibrium SMS is uniformly 
rotating initially, it cannot support a large amount of 
angular momentum without exceeding the mass-shedding 
limit. However, when internal magnetic fields and viscos- 
ity are weak, the equilibrium star may rotate differentially 
and support a considerably larger angular momentum (and 
higher T/W) without exceeding the mass-shedding limit. 
For such a case, New & Shapiro (2001) propose an alter- 
native scenario for the quasi-static cooling and contraction 
of an equilibrium SMS which inevitably leads to bar for- 
mation prior to catastrophic collapse. Such a scenario will 
almost certainly generate long wavelength quasi-periodic 
wave emission, even in the absence of SMBH formation. 

Our (3+1) hybrid PN calculations offer a first glance at 
SMS collapse. An improved description will require sev- 
eral refinements to our computational scheme. First, it 
will be necessary to employ a fully relativistic treatment 



of Einstein's field equations to explore the final dynamical 
phase of collapse once a BH has formed. In fact, we can- 
not reliably use PN gravity once the radius decreases much 
below Rp/M < 20, since the central fields are becoming 
strong. One possibility would be to use our output when 
the central fields are beginning to get strong (a c < 0.5) as 
initial data for a fully relativistic (3+1) code. This would 
allow us to take advantage of the faster, more robust PN 
code to handle the Newtonian and PN implosion regime 
and switch over to a fully general relativistic code only for 
the final strong-field phase of collapse. 

Secondly, our computation would benefit significantly 
from nested grid method (see, e.g., Ruffert 1992) to han- 
dle the large dynamic range characterizing SMS collapse. 
We have successfully exploited approximate homology by 
introducing a scale factor to set up an approximate comov- 
ing coordinate system. But homology breaks down at the 
center toward the end of the collapse. Moreover it is diffi- 
cult to follow the small amount of matter in the equatorial 
plane that is supported by centrifugal forces while the bulk 
of the matter collapses. Nested grid method may provide 
one means of concentrating computational resources on 
the central region of the star while simultaneously resolv- 
ing the low density outermost regions. 
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